Force balance in canonical ensembles of static granular packings 
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We investigate the role of local force balance in the transition from a microcanonical ensemble of 
static granular packings, characterized by an invariant stress, to a canonical ensemble. Packings in 
two dimensions admit a reciprocal tiling, and a collective effect of force balance is that the area of this 
tiling is also invariant in a microcanonical ensemble. We present analytical relations between stress, 
tiling area and tiling area fluctuations, and show that a canonical ensemble can be characterized by 
an intensive thermodynamic parameter conjugate to one or the other. We test the equivalence of 
different ensembles through the first canonical simulations of the force network ensemble, a model 
system. 



PACS numbers: 45.70.Cc, 05.40.a, 46.65. +g 

Dense packings of grains are athermal: having achieved 
a mechanically stable state, they remain there unless ex- 
ternally driven. The set of final states of a particular 
numerical or experimental preparation protocol defines 
a noncquilibrium ensemble of static packings. A funda- 
mental and open question, ultimately related to the pack- 
ing and mechanical properties of granular media, is the 
proper framework for the statistical description of such 
an ensemble. This work addresses one important aspect 
of the problem, namely the role of local force balance 
and how it influences the transition from a microcanoni- 
cal ensemble, here defined as an ensemble at fixed stress, 
to a canonical ensemble. 

Energy is a natural quantity to characterize equilib- 
rium states because it is an invariant of the dynamics. 
This distinction is lost in an ensemble of static pack- 
ings, and one must search for other convenient quantities 
to replace it. A number of prior works have proposed 
characterizing a statistically homogeneous packing by its 
average stress [H O El HI El El O [9] . This is a natural 
choice in the sense that the "extensive stress" aV [3D] on 
a body of volume V in mechanical equilibrium under an 
external load is determined by the forces on its boundary, 
regardless of the configuration of grains in the bulk |10j . 
A microcanonical ensemble at fixed extensive stress may 
then be defined as the set of all arrangements of N grains 
consistent with a particular boundary loading. Through- 
out this work we restrict ourselves to frictionless packings 
of disks in two dimensions and to isotropic stress states, 
so that the extensive pressure V = ^(Trer)P suffices to 
characterize the stress. The extensive pressure is addi- 
tive: V — ^2ipi- In a disk packing the "pressure" on a 

grain is pi = ^ • fij ■ fij , where fij , the force that grain 
j applies to grain i, is nonzero only if the grains are in 
contact and fij is the vector from the center of j to i. 

Postulating entropy maximization, a Boltzmann-likc 
factor exp (—aV) follows for the canonical ensemble, pro- 
vided V is the only relevant invariant of the microcanon- 
ical ensemble. The quantity a is an intensive thermody- 
namic parameter conjugate to V . Edwards has recently 



suggested the name angoricity for a~ l [S]. It is distinct 
from the more well-known compactivity [11 , which is 
conjugate to the packing volume V. 

To date, the main application of stress-based ensembles 
of static packings has been the prediction of the statis- 
tics of local measures of stress, such as the force / at 
a contact or the pressure p on a grain, which provide a 
fundamental characterization of the stresses in a pack- 
ing. The forces on one grain are coupled to the forces 
on other grains via Newton's third law. In a statistical 
treatment, local force balance enters through (^-functions 
in the partition function of the system |12j , which must 
then be integrated over; in analytical calculations this is 
tedious or impossible for systems larger than a few grains 
[121 fT3] . All published calculations resort to some degree 
of approximation to treat local force balance, and most 
neglect spatial correlations completely. At the latter level 
of approximation, the tail of a local stress probability 
distribution function reflects the form of the Boltzmann 
factor, and hence these approaches predict local stress 
probability distribution functions that decay exponen- 
tially for large stresses [TJ HI El [14] . This prediction 
is analogous to the Maxwell-Boltzmann distribution in 
an ideal gas, which is exponential in the particle energy. 
We will refer to calculations that neglect spatial corre- 
lations as "ideal gas- like" , though this is not meant to 
invoke a gaseous state. 

In recent work [15j . we pointed out that it is possible 
to improve on ideal gas-like calculations in 2D packings 
by making use of a dual structure known as the Maxwell- 
Cremona diagram or reciprocal tiling [16] . In Section |TJ 
we construct the reciprocal tiling and explain how it in- 
fluences local stress statistics. Crucially, the tiling exists 
as a necessary consequence of local force balance. The 
most important feature of these tilings is their area A, 
which we shall see is an extensive invariant much like V '. 
As the constraint of mechanical equilibrium is one of the 
principal differences between ensembles of granular pack- 
ings and other ensembles, a key question is how the tiling 
area A can be incorporated in an ensemble treatment of 
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FIG. 1: (a) A network of balanced contact forces on the 
frictionless periodic triangular lattice. Nodes correspond to 
grains and edges to contacts. Color and line thickness are 
mapped to contact force magnitude. The contact network 
has dimensions L\ x L2. The total force on a surface parallel 
to the boundary of the unit cell with dimension L\ (L2) is F\ 
(F2). (b) A single-grain state in the triangular lattice. Ar- 
rows represent vector contact forces acting on the grain, (c) 
The reciprocal tiling, or Maxwell-Cremona diagram, corre- 
sponding to the force network in (a). Each tile is constructed 
from the contact forces acting on a grain. The unit cell of the 
tiling has dimensions F\ x F2. (d) Construction of a tile for 
the grain in (b). Vector forces, rotated by 7r/2, are graphically 
summed around the grain; the tile closes because the grain is 
in static force balance. The tiles tessellate space due to New- 
ton's third law. Rotating forces by -n/2 is not essential, but 
doing so ensures that, e.g., if grain j is to the right of grain i 
in the packing, tile j will be to the right of tile i in the tiling. 



static granular packings beyond the ideal gas approxi- 
mation. The remainder of this work seeks to answer this 
question. 



I. RECIPROCAL TILINGS 

The Maxwell-Cremona diagram or reciprocal tiling is 
a dual structure constructed from the forces in a two- 
dimensional packing; it is a geometric representation of 
the stress state (see Fig.[l^i,c). Each grain maps to an in- 
dividual tile in the tiling. The boundaries of this tile are 
constructed by graphically summing the contact forces 
on the grain, moving from contact to contact around the 
grain in a right hand fashion (see Fig. [lj>,d). A tile's 
boundary closes because the vector sum of the forces on 
the grain is zero, i.e. because the grain is in force bal- 
ance. Moreover, tiles tessellate space due to Newton's 
third law, which guarantees tiles in contact have facets 
with like length and orientation. Note that lengths in the 
tiling correspond to forces in the packing, so the tiling oc- 
cupies a different space than the packing. In particular, 
the area of the tile corresponding to grain i has units 
of (force) 2 . Modulo a global rotation, the vertex coordi- 



nates in a Maxwell-Cremona diagram are equivalent to 
the "void forces" of Satake [17] or "loop forces" of Ball 
and Blumcnfeld [15] , 

As noted above, the set of grain arrangements com- 
patible with fixed boundary load specifies an ensemble 
at fixed extensive stress aV [7]. Similarly, the boundary 
of a packing's tiling can be constructed simply by know- 
ing the boundary forces. If we then rearrange the grains 
or forces inside the packing to produce a new static pack- 
ing compatible with the same boundary loading, the the 
tilings corresponding to the old and new packings will 
have the same area. Therefore the total area A = J2i a i 
of the Maxwell-Cremona tilings in an ensemble at fixed 
extensive stress is also an additive invariant |15l 119] . The 
invariance of A is a direct consequence of local force bal- 
ance in the packing and has important consequences for 
the statistics of local stresses. 

In Rcf. 15] we considered an ideal gas-like calcula- 
tion in which the average tiling area (A), as well as the 
average extensive pressure (V), is imposed. Maximiz- 
ing entropy then leads to a local pressure distribution 
P(p) = Z~ x p v exp (— ap — 7(a(p))), where v depends on 
contact geometry and Z, a, and 7 are Lagrange multipli- 
ers. The distribution is asymptotically Gaussian because 
(a(p)), the average tile area given p, is proportional to 
p 2 . As shown in Fig. [2J this distribution is in excellent 
agreement with numerics in a model system, the force 
network ensemble of Snoeijer et al. [20], and a clear im- 
provement over an ideal gas-like calculation that does 
not enforce (A). As we discuss further below, the suc- 
cess of this approach is suggestive of a Boltzmann factor 
exp (— aT — jA) that incorporates A in addition to V . 

The form of the tail of local stress distributions in 
static granular packings is a subject of ongoing debate, 
and our goal here is not to insist on one form over an- 
other. Rather, we ask what local stress distributions can 
tell us about applying the maximum entropy postulate 
to ensembles of packings. Here the force network ensem- 
ble is particularly useful, as it serves as a litmus test for 
theory. A calculation that is too simplistic to explain 
results in the force network ensemble cannot explain re- 
sults in numerical or experimental ensembles, which are 
more complex than the force network ensemble; appar- 
ent agreement, if any, must be coincidental. Therefore, 
ideal gas-like calculations that predict exponential tails 
cannot explain the exponential tails observed in some ex- 
perimental and numerical measurements [211122] . because 
the same calculations should apply to the force network 
ensemble, which does not display exponential tails. 

As the tiling area constraint is a necessary conse- 
quence of force balance and incorporating the constraint 
in an ideal gas-like calculation yields predictions consis- 
tent with the force network ensemble, it is important to 
ask how the tiling area should enter, more generally, in a 
maximum entropy approach. A more detailed calculation 
than that in Ref. [15] is likely required to describe sys- 
tems more complex than the force network ensemble, for 
example due to a growing correlation length 23 , 24J . The 
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FIG. 2: Numerical and theoretical probability distribution 
functions (see legend) of local pressure P(p) in the force net- 
work ensemble on a frictionless triangular lattice; adapted 
from Ref. |15| . The numerical distribution is taken from a 
microcanonical ensemble of N = 1840 grains using umbrella 
sampling [251 126] . Theoretical distributions result from an 
entropy maximization calculation subject to a constraint on 
the average extensive pressure (V) and/or average area of a 
reciprocal tile (A). All calculations neglect correlations with 
neighboring grains, while the numerics impose force balance 
exactly on every grain. 



results of Ref. [15], then, suggest two possible scenarios, 
but do not go far enough to distinguish between them. 
These are: (i) The extensive quantities V and A must be 
treated on equal footing, i.e. the Boltzmann factor should 
in fact be exp (— aV — jA). This possibility was already 
noted above. The second is that (ii) the tiling area A 
need not be enforced independently, that is it need not 
appear in the Boltzmann factor. Instead, its role is that 
of proxy for the spatial couplings implied by local force 
balance and neglected in an ideal gas calculation, and if 
these are incorporated exactly in an analytical calcula- 
tion the tiling area A need not be considered separately. 
The goal of the present work is to distinguish between 
these two scenarios, which is a prerequisite for any future 
work that would seek to incorporate spatial correlations 
in more detail. The issue hinges fundamentally on the 
role of local force balance in a static granular ensemble. 

By considering two routes by which a microcanonical 
ensemble passes to a canonical ensemble, we argue that V 
and A are not independent in the thermodynamic limit, 
i.e. we argue in favor of the latter of the two scenarios 
above. As support we offer the first analysis of the canon- 
ical force network ensemble. We will show that for small 
systems it matters a great deal which Boltzmann factor 
is used. For a thermodynamically large system, however, 
it suffices to characterize the system by either V or A, 
provided one also imposes force balance on every grain. 
We confirm this through simulations of the force network 
ensemble, which allows us to access the canonical ensem- 
ble directly with Monte Carlo methods and to impose 
local force balance exactly. 



II. FORCE NETWORK ENSEMBLE 

The force network ensemble [2D] is an ideal testbed for 
statistical mechanics-based approaches to static granu- 
lar media. We frame our discussion in the context of 
the force network ensemble, but our main conclusions re- 
garding ensembles based on V and/or A apply to any 
ensemble of static granular packings. For now it suf- 
fices to state that, rather than comprising many different 
arrangements of the grains, the force network ensemble 
takes advantage of the fact that disk packings at finite 
pressure are generically hyperstatic, i.e. the forces are un- 
derdetermined by the constraints of mechanical equilib- 
rium. The ensemble then comprises all configurations 
of forces (force networks) on one quenched configuration 
of grains. In the microcanonical force network ensem- 
ble the global stress tensor is imposed, and all force bal- 
anced configurations of noncohesive forces are assigned 
equal statistical weight. For purposes of illustration we 
take the contact network to be a frictionless triangular 
lattice of grains and discuss at several points the relation 
to disordered packings. The force network ensemble is 
described in greater detail in Section |lVj 



III. PASSING TO THE CANONICAL 
ENSEMBLE 

As a thought experiment, one passes from a micro- 
canonical to a canonical ensemble by placing a previously 
microcanonical system in contact with a much larger sys- 
tem that acts as a reservoir of some conserved quantity, 
e.g. V. For ensembles of static packings, this can be 
achieved by sampling A-grain clusters of grains from a 
microcanonical ensemble of packings of M 3> N grains. 
This is illustrated in Fig. (3]i and b. The extensive pres- 
sure in the canonical system V v fluctuates from clus- 
ter to cluster, but its average is dictated by the bath: 
(V u ) = N(p), where (p) is the average pressure per grain 
in the bath. Formally, the parameter a in the Boltzmann 
factor exp (— aP v ) plays the role of a Lagrange multiplier 
that imposes this constraint on the canonical system. 

The twist is that, as noted above, a microcanonical sys- 
tem at fixed V also has fixed tiling area A. Therefore the 
bath also imposes an average tiling area (A v ) — N(a) on 
the canonical system, where (a) is the average area per 
tile in the bath. If the imposed average tiling area is inde- 
pendent of the imposed average extensive pressure, then 
an additional Lagrange multiplier 7 is needed, i.e. the 
Boltzmann factor should be exp — jA u )- 

To illustrate the potential complications introduced 
by A in a canonical ensemble, consider the non-periodic 
canonical system depicted in Fig. [1] (We now dispense 
with the subscript v for extensive quantities in canonical 
systems.) The system with dimensions Lx L experiences 
a loading due to forces imposed by the bath. In terms 
of its boundary forces {/ c } acting at positions {x c }, the 
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FIG. 3: (a) A periodic triangular lattice. In the force net- 
work ensemble, the set of force networks (see Fig. [T^,) with 
fixed stress tensor comprise the microcanonical force network 
ensemble on this fixed contact network, (b) A canonical en- 
semble can be created by embedding the system of (a) in a 
larger packing that then acts as a bath. The extensive stress of 
the system and bath, together, is conserved, and the canoni- 
cal system is no longer periodic, (c) Alternatively, the system 
of (a) can be placed in contact with a bath, i.e. allowed to 
exchange extensive stress with the bath in such a way that 
the extensive stress of the two systems is conserved and the 
canonical system remains periodic. 

system's extensive stress is [TU] [3T] 

CafiV = f c ,a X Ct/ 3 . (1) 

c 

For boundary forces of magnitude / at orientation 9, 
separated by a distance d as in Fig. [i] the system's 
extensive stress tensor is isotropic, with extensive pres- 
sure V = 2 f(L cos 9 + dsm9). The tiling area is A — 
2/ 2 (l + sin 29 + cos26>). Note that A can be changed 
while holding V fixed by varying / and 9 and requiring 
d= {V-2fL cos 6)/2f sintf. 

The above example suffices to show that A cannot be 
a single-valued function of V in a non-periodic system. 
Introducing more, and more disordered, boundary forces 
will only enhance this degeneracy, as there will be many 
ways to choose the boundary forces to achieve a particu- 
lar V, and these configurations will have different tiling 
areas. The question is whether the relative area fluctu- 
ations W {8 A 2 ) I (A) at fixed V become negligible in the 
thermodynamic limit; if so, imposing (V) via the inten- 
sive parameter a also suffices to select (^4). If not, an ad- 
ditional thermodynamic parameter conjugate to A must 
be introduced. 



A. Periodic canonical systems 

Before returning to the fluctuations of tiling area in 
the system of Fig.|3]D, it is useful to consider an alternate 
route from the microcanonical to the canonical ensemble. 
While in Fig. :3jb the canonical system is non-periodic and 
embedded in a larger system, we now imagine placing a 
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FIG. 4: (a) Boundary loading of a square body, (b) The 
boundary of the reciprocal tiling of the system in (a). With 
the appropriate variation of boundary force magnitude /, ori- 
entation 6, and separation d, the area of the tiling can be 
changed while the stress is held fixed. 

periodic system in contact with a reservoir with which it 
is allowed to exchange extensive stress, as illustrated in 
Fig. [3]:. Unlike the previous scenario, which can be un- 
derstood as clusters sampled from a larger microcanoni- 
cal ensemble, this scenario corresponds to a thought ex- 
periment in which a collection of N grains is repeatedly 
randomly seeded inside a prescribed unit cell and allowed 
to relax to mechanical equilibrium. Final states that 
are mechanically stable are placed in the ensemble and 
weighted according to the Boltzmann factor. 

The key observation is that, in a periodic force net- 
work, the tiling area is a single- valued function of the ex- 
tensive stress tensor: A = (deta)V. To see this consider 
the force network and tiling in Fig. [T]i and c. Imposing 
the stress tensor is equivalent to imposing the net vector 
forces F\ and F2 acting on two boundaries of the L\ x L2 
unit cell. For the force network in Fig. [T]i the extensive 
stress tensor is 

where V — L\L2- It is straightforward to see that, due 
to periodicity, the area of the network's reciprocal tiling 
must be A = FiF 2 = (det<r)V". As one can always con- 
struct a rectangular unit cell and choose coordinate axes 
aligned with the principal stress directions, the relation 
A = (deter) V is general for periodic force balanced net- 
works in two dimensions, independent of contact network 
geometry and topology. 

The relation between tiling area and stress has imme- 
diate consequences. First, in a canonical ensemble of 
periodic packings (Fig. [3J;) it is clear that extensive pres- 
sure V and tiling area A are not independent; knowl- 
edge of V implies knowledge of A. For ensembles of 
noncohesive isotropic packings, the relation is bidirec- 
tional: V can be inferred from A. Therefore it suffices to 
employ a Boltzmann factor exp {—aP) or, alternatively, 
exp (— 7-4). This latter possibility has not previously 
been considered, and we test it below. 

The second consequence follows from noting that the 
routes to the canonical ensemble depicted in Fig. and 



c differ only in their implementation of boundary con- 
ditions. In general, one anticipates that details at the 
boundary should not influence the thermodynamic limit, 
and thus in the thermodynamic limit a Boltzmann fac- 
tor exp {—aP) or exp (—jA) will also suffice for non- 
periodic systems. Consequently, it is not necessary to 
employ a Boltzmann factor exp (— aP — jA) in a thermo- 
dynamically large system. We reinforce this expectation 
below by providing a scaling argument for the relative 
area fluctuations in a non-periodic system. We will test 
these predictions by performing simulations of the canon- 
ical force network ensemble for both periodic and non- 
periodic force networks of varying size N. We emphasize 
that equivalence in the thermodynamic limit does not 
imply equivalence in small systems, a point we discuss 
further below. 



B. Non-periodic canonical systems 

Returning to non-periodic canonical systems, as in 
Fig. [HJa, we consider the tiling area fluctuations in the 
thermodynamic limit. We argue that for each extensive 
pressure V the system approaches the same tiling area 
selected in periodic systems with vanishing relative fluc- 
tuations. 

Consider a square system subject to isotropic compres- 
sive loading. One of its boundaries is depicted in Fig. [5^i. 
The boundary is subject to a net force (0,—F), with 
O(VJV). The iV b - 0(VN) boundary forces 
can be used to construct one boundary of the system's 
reciprocal tiling, shown in Fig. [5Jd. This boundary re- 
sembles a directed random walk starting from the origin 
and constrained to end at (F, 0) . Each possible walk 
is as likely as its reflection about the dashed segment 
in Fig. [51). It follows that the mean tiling area will be 
(A) = F~ = (det<r)V, just as in a periodic system. 

It remains to consider the dependence on N of the 
tiling area fluctuations •*/ {5 A 2 ). We assume each step 
of the walk goes a typical distance (/) = F/N^ to the 
right in Fig. [5J). In the following the load per unit length 
(/) = F/Nh is kept constant as the system size increases. 
Labeling hi the y-coordinate of the walker at step i, the 
typical vertical distance traveled after i steps is \J (h 2 ). 
Sufficiently close to the left endpoint of the walk, the 
influence of the constraint to end at (F,0) will not be 
felt, and we must have \J (h 2 ) ~ (f)Vi- Similarly, near 
the right endpoint \J (h 2 ) ~ (f)y/Nb — i- 

We want to develop an upper bound on the tiling area 
fluctuations. In this spirit, we assume that the typical 
vertical displacement behaves like a simple random walk 
up to the midpoint, i.e. 
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(f)Vi i < lN h 

(f)VW^i i>lN h . 
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The scaling of absolute area fluctuations under a simple 
random walk is captured by J (S A 2 ) ~ (/) J2i V C 1 *)- 



FIG. 5: (a) Boundary of a non-periodic structure subject to 
externally imposed forces. The net force (-F, 0) is purely 
compressive, (b) The boundary of the structures reciprocal 
tiling can be imposed from its boundary forces. It resembles a 
directed random walk constrained to travel a distance (F, 0). 
Its vertical displacement after i steps is labeled hi. 

Using Q we find y^JjA 2 ) < </) 2 iv| ~ O(iVf). As 
(A) is extensive, the relative fluctuations \J (5A 2 ) /(A) 
decay at least as fast as N~i, and therefore vanish in 
the thermodynamic limit. 

Note that our argument relies on treating boundary 
forces as random variables. For a packing at isostatic- 
ity this cannot hold; specifying only half the boundary 
forces of an isostatic packing suffices to fix the other half 
[27] , We now show that if a packing has mean contact 
number z = Zj SO + Az, there is a length scale £* such 
that, for packings of linear dimension L > £*, all the 
boundary forces can indeed be treated as random vari- 
ables and the above scaling argument can again be in- 
voked. In order to randomly assign all the JVb boundary 
forces, there must be enough underdetermincd forces in 
the bulk. There are N c ~ Az N of these "excess" forces. 
For N e to balance or exceed iVj, requires Az> AT- 5. As 
the system's linear dimension L ~ %/iV, this is equivalent 
to requiring L > l/Az. Note that the same balance be- 
tween boundary contacts and excess bulk contacts is also 
used to derive the isostatic length £* ~ l/Az [23l [24] . 
which governs the crossover from discrete to continuum 
response in static packings. Therefore our argument pre- 
dicts that relative area fluctuations vanish in the thermo- 
dynamic limit for systems arbitrarily close to isostaticity, 
provided the system size L exceeds the diverging length 
scale £* ~ l/Az. Our scaling argument cannot be ap- 
plied to systems smaller than £*, so it may be possible to 
identify stronger bounds on the area fluctuations. 



IV. NUMERICAL RESULTS 

A. Monte Carlo methods 

The force network ensemble is a convenient venue to 
test the ideas described above. Because the ensemble can 
be sampled with Monte Carlo methods, it is possible to 
realize, numerically, the thought experiments described 
in Fig. [3] 

We briefly summarize the force network ensemble and 
the Monte Carlo methods we employ to sample the 
canonical ensemble. Having selected a hyperstatic con- 
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tact network, the key requirement is the imposition of 
force balance on each grain i, i.e. J2j fij = 0' where 

fij is the contact force on i due to j. The micro- 
canonical force network ensemble also imposes the stress 
ct q /3 = (r/V)J2tj hj n ij,anij,p, where r is the grain ra- 
dius, V is the area of the quenched packing, and n, 3 is the 
unit vector from i to j. We restrict ourselves to isotropic 
systems, and the important quantity will be the extensive 
pressure V = (V/2)Tr<7 = (r/2) J2ij fij- We also impose 
a positivity condition on each force, /y > 0, restricting 
the ensemble to noncohesive force networks. Any config- 
uration of forces f = {fij} satisfying these constraints is 
an element of the microcanonical force network ensemble; 
Fig. [T]i gives an example. 

The ensemble can be efficiently and flatly sampled us- 
ing Monte Carlo techniques. For periodic systems we use 
the "wheel moves" of Ref. [T3] . These may be modified to 
work in non-periodic systems, as well [25, 26J. Changes to 
boundary forces are allowed provided they respect force 
balance and noncohesivcncss of the forces. In the canon- 
ical ensemble we add an additional move that changes 
the extensive pressure V . In a triangular lattice, this 
move simply adds a quantity e to each force. A proposed 
move from a force network f„ to a force network f v i is 
accepted with a probability determined according to the 
Metropolis acceptance rule, 

acc(f„ -> f„) = min (l, 0(f„) . (4) 

Here B(f v ) is the Boltzmann factor, which we vary below. 
0(f„') = 1 if the configuration f„/ is force balanced and 
noncohesive and otherwise. 



B. Microcanonical ensemble 

All previous simulations of the force network ensemble 
have been of a microcanonical variety [HI [T31 E3 HH1 [201 
|25j|26]. The microcanonical frictionless triangular lattice 
has been well-studied in the force network ensemble, and 
local stress statistics have been determined extremely ac- 
curately by employing umbrella sampling |15L 125] . Fig. [6] 
plots the microcanonical P(p) for systems of N = 115, 
460, and 1840 grains, demonstrating that finite size ef- 
fects are negligible in the largest system. 

C. Canonical ensembles 

To test the above predictions we perform simula- 
tions of the canonical force network ensemble on peri- 
odic and non-periodic triangular lattices. We consider 
two ensembles equivalent if their local pressure distri- 
butions P(p) converge to the same limiting distribution 
as N — ► oo. We take advantage of the fact that P(p) 
is known extremely accurately in the microcanonical en- 
semble (Fig.pl and consider P(p) in a microcanonical cn- 




FIG. 6: Local pressure probability distribution P(p) in the 
microcanonical force network ensemble on a periodic friction- 
less triangular lattice for unit cells containing N = 115, 460, 
and 1840 grains. Finite size effects are negligible. 
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FIG. 7: (a) Local pressure probability distribution P(p) in 
the canonical force network ensemble in periodic systems with 
Boltzmann factor exp (—otP) (solid curves) and a microcanon- 
ical system with N = 1840 (dashed curve). The canonical 
system size N = 18, 24, 30, 56, 120, and 224 increases in the 
direction of the arrow, (inset) The parameter a is insensi- 
tive to system size and consistent with a{p) = 1. (b) P(p) in 
the canonical ensemble in periodic systems with Boltzmann 
factor exp (—7.4) (solid curves) and microcanonical ensemble 
(dashed curve). System sizes are identical to (a), (inset) The 
parameter 7 is insensitive to system size and consistent with 
7(a) = |. 



semble of N — 1840 grains representative of the thermo- 
dynamic limit. By studying canonical systems of varying 
system size N, we look for numerical evidence of con- 
vergence to P(p) from the large microcanonical system. 
By checking if different canonical ensembles converge to 
this distribution, we can compare them to each other by 
transitivity. We emphasize the perspective that a and 
7 are Lagrange multipliers: for systems of different size 
N, thermodynamic parameters are not assigned but var- 
ied in order to achieve a predetermined {V)/N and/or 
(A)/N. 

We first test the prediction that a periodic canonical 
ensemble can be described using either of two Boltzmann 
factors. Fig. [7^ depicts the case B(f v ) = cxp(— aV v ). 
Systems of varying size N are simulated, and a is cho- 
sen so that (V) — 6N. The probability distribution P(p) 
of the local pressure p is sampled, and converges to the 
distribution sampled from a large microcanonical force 
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network ensemble, a is nearly constant over the range of 
system sizes sampled and consistent with its expected 
value in the thermodynamic limit a(p) — 1 (see Ap- 
pendix) . This strongly indicates that the microcanonical 
force network ensemble and the canonical force network 
ensemble with Boltzmann factor exp (— aP) are equiva- 
lent in the thermodynamic limit. It is remarkable that 
for small systems the asymptotic decay of P(p) appears 
nearly exponential, though this is not the behavior in 
the thermodynamic limit. As always, caution is required 
when interpreting simulations or measurements in small 
systems. 

Fig. [rj^ depicts the case B({ u ) = exp(— jA u ) applied 
to an ensemble of periodic force networks. The proba- 
bility distribution P(p) is again sampled for varying sys- 
tem size. The intensive parameter 7 is chosen to en- 
sure {A} = ^-N, which follows from A = (deter) V for 
(V) = 6N. The pressure statistics again converge to the 
microcanonical distribution in Fig. [6] and the thermody- 
namic parameter 7 is independent of N and consistent 
with its expected value 7(a) = | (see Appendix). This 
provides strong evidence that the microcanonical force 
network ensemble and the canonical force network en- 
semble with Boltzmann factor exp (—"/A) are equivalent 
in the thermodynamic limit. By transitivity it also indi- 
cates equivalence of the two canonical ensembles in the 
thermodynamic limit. 

We next consider non-periodic systems. Packings are 
composed of grains organized in hexagonal layers around 
a central grain. P(p) is sampled on the central grain. 
We employ a Boltzmann factor B{f u ) = exp (— aV — ^A) 
and choose a and 7 such that (V) = 6A and (A) = 
^Y^N. If the two constraints are not independent, one 
of the thermodynamic parameters must tend to zero. 
The thermodynamic parameters a and 7 are plotted in 
Fig. [8J). Indeed, 7 tends to zero with increasing system 
size, consistent with argument that relative tiling area 
fluctuations vanish with increasing N. 

Though a in Fig. exceeds its value from periodic 
packings (Fig. [7^ inset), we argue that it will ultimately 
approach the same limiting value. It may be shown (see 
Appendix) that the two parameters satisfy the relation 



a(p) + 2 7 (a) 



Nc 
N 



(5) 



where N c is the number of contact forces in the sys- 
tem. Contours of (|5| are plotted in Fig. [8] The ap- 
proximations implicit in (|5| become increasingly accu- 
rate as the system size increases; indeed, the largest sys- 
tem sizes in Fig. [8] lie on their respective contours. The 
intercept value — 2 of the contours provides an up- 
per bound on a(p) and smoothly approaches its limiting 
value j± — 2 — ► \{z — Z{ so ) from above, where z — 6 
is the mean contact number and Zj SO is the isostatic co- 
ordination number for frictionless disk packings. Hence, 
although a(p) overshoots its asymptotic limit, we antici- 
pate that it will approach the asymptotic value smoothly 
for system sizes larger than we can simulate practically. 




1 a(p) 



FIG. 8: (a) Local pressure probability distribution P(p) in 
the canonical force network ensemble for non-periodic sys- 
tems with Boltzmann factor exp (—aP — "/A) (solid curves) 
and a periodic microcanonical system (dashed curve). The 
canonical system size (N — 7, 19, 37, 61, and 91, correspond- 
ing to 1, 2, 3, 4, and 5 hexagonal layers around a central grain, 
respectively) increases in the direction of the arrow. The mi- 
crocanonical distribution is identical to that in Figs. [2] and 
[7] (b) Evolution of the parameters a and 7 (triangles) with 
system size. Lines are contours of |5| for the same systems 
sizes, as well as the thermodynamic limit. 



It is evident from Figs. [7] and [8^1 that local pressure 
statistics in small systems are closer to their asymp- 
totic form in the periodic system with Boltzmann factor 
exp (—7^4) than in the periodic system with Boltzmann 
factor exp (— o/P). For the smallest systems they are com- 
paratively closest to their asymptotic form in the non- 
periodic system with Boltzmann factor exp (— aV — "/A). 
This is quantified in Fig. [9j which plots the convergence 
of local pressure fluctuations A N := \J {Sp 2 ) n / {p) N to 
their N — ► 00 value, approximated by Ajv evaluated for 
the numerically sampled microcanonical distribution for 
N = 1840. For all cases, A^ decays with N, confirm- 
ing the convergence of canonical and microcanonical en- 
sembles that was already apparent from Figs. [7] and [8] 
For the system sizes sampled in both periodic and non- 
periodic systems, Ajv is always smaller when the inten- 
sive thermodynamic parameter 7 is employed. In peri- 
odic systems (A at — ^oc) — A -1 , while in non-periodic 
systems the decay with system size is weaker. 



V. CONCLUSIONS 

Static packings possess reciprocal tilings as a conse- 
quence of local force balance. The area of the tiling is 
therefore an additive extensive quantity that is induced 
by the presence of local force balance. We have con- 
sidered the role of this tiling area, and through it force 
balance, in the construction of a canonical ensemble of 
static granular packings. 

In a canonical ensemble both the extensive stress and 
the tiling area fluctuate from configuration to configu- 
ration. These fluctuations are not independent: view- 
ing the bath as a microcanonical system in which the 
canonical (sub)system is embedded, the bath pressure 
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FIG. 9: Evolution of the relative fluctuations Ajvwith sys- 
tem size. Fluctuations from finite size systems are taken from 
(a) the periodic systems of Fig. u\ with Boltzmann factors 
exp (— aP) and exp (—7.4), and (b) the non-periodic systems 
of Fig. [8] with Boltzmann factor exp (— aV — 7-4), as well 
as non-periodic systems sampled with the Boltzmann factor 
exp (— aP) (distributions of P(p) not shown). The fluctu- 
ations in the thermodynamic limit are estimated from the 
distribution in Fig. [2] 



V = ^(Tra)V and tiling area A = (deta)V are related. 
Imposing an average extensive pressure (V) in the canon- 
ical system, via a conjugate intensive parameter a, also 
suffices to select the proper average tiling area (A) in the 
thermodynamic limit. This expectation was confirmed 
both with a scaling argument for the relative tiling area 
fluctuations and numerical simulation. 

Note that our simulations directly addressed only or- 
dered packings far from isostaticity. Recalling that the 
relation between A and V in periodic systems holds re- 
gardless of disorder, there is no reason to expect disorder 
to prevent convergence in the thermodynamic limit. A 
disordered system must of course be large enough to aver- 
age over fluctuations in the local contact number; there- 
fore, ceteris paribus, a numerical measure of convergence 
such as Ajv will have a smaller prefactor for ensembles 
of ordered packings than for disordered packings. We ar- 
gued above that in the limit Az J, 0, convergence to the 
thermodynamic limit can only be expected for systems 
of size L ^S> i* ~ 1/Az. This is consistent with previous 
studies [21] but should be tested explicitly. Numerical 
tools beyond the force network ensemble are required, 
however, because the ensemble itself vanishes in the iso- 
static limit. 

We have also demonstrated that a Boltzmann factor 
exp (— 7-A) may be used in place of exp {—aP). This may 
seem a curiosity; e.g. the equilibrium Boltzmann factor 
exp (—(3E) can also be replaced with exp (—f3E 2 ) in a 
purely repulsive system, though additivity is lost and no 
insight is gained. There are important differences, how- 
ever, for ensembles of static packings. Unlike E 2 , the 
tiling area is an additive quantity, just as the extensive 
pressure. More importantly, the tiling area adds new 
insight by relating mechanical equilibrium to geometry. 
As an additional practical consideration, stress statistics 
converge faster when the intensive thermodynamic pa- 
rameter is conjugate to A. Finally, a Boltzmann factor 



of the form exp (—7.4) supplies the physical intuition, 
confirmed in the force network ensemble, that local stress 
statistics in two dimensions should display Gaussian tails. 
In a system with Boltzmann factor exp (— aP), an ideal 
gas-like approximation incorrectly predicts exponential 
tails. We have demonstrated numerically that P(p) in a 
canonical system with Boltzmann factor exp (—aP) and 
local force balance does indeed approach the correct form 
with a Gaussian tail for sufficiently large system size. To 
reproduce the Gaussian tail in a calculation employing 
this Boltzmann factor, however, would require integrat- 
ing over local force balance constraints in an "interact- 
ing" (via Newton's third law) system [25] , 

In analytical calculations, it is often difficult or im- 
possible to integrate out local force balance constraints 
in a large, potentially disordered system. Approxima- 
tions to local force balance, e.g. neglecting spatial corre- 
lations or using scalar force balance |14j , break the mech- 
anism whereby imposing (V) also produces the correct 
(.4). In this situation, it becomes helpful to simultane- 
ously impose both pressure and tiling area constraints 
in an entropy-maximization calculation, as in Ref. |15j . 
In effect one demands that the system satisfy two global 
constraints, one of which would have followed "for free" 
due to local constraints, i.e. force balance, had they been 
incorporated exactly. Parameters a and 7 determined 
this way can no longer be identified with their values 
in the thermodynamic limit; instead they should simply 
be viewed as Lagrange multipliers imposing a constraint. 
The resulting gain is a dramatic improvement in the ac- 
curacy of the predicted statistics, as in Fig. [2] even in an 
ideal gas-like calculation. 
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APPENDIX A: INTENSIVE THERMODYNAMIC 
PARAMETERS 

Assuming a Boltzmann factor exp (— aP — 7^4), we can 
calculate a relationship in the thermodynamic limit be- 
tween the thermodynamic parameters a and 7. We take 
the system to have N grains and N c contacts. The par- 
tition function is 

Z = J die- aV ^ A ^Q(i) 

dPdAG(P,A) e- aV - lA . 



(Al) 



6(f) restricts to balanced noncohesive force networks. 
G(P, A) is the density of states with extensive pressure 
V and tiling area A, which can be rewritten G(P,A) = 
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D.(V)^(A\V). nCP) is the density of states with V, which 
scales as fi('P) ~ pN c -2N-i m a f r j c tionless disk pack- 
ing [T5]. ^/(^IT 5 ) is the conditional density of states with 
tiling area A given extensive pressure V. The thermody- 
namic parameters a and 7 must be chosen to ensure 



(V) 
(A) 



dlnZ 

da 
dlnZ 



(A2) 



The conditional probability W is sharply peaked near 
(A^)), the average area of a tiling with extensive pres- 
sure V jlH]- As the system size grows, this function is in- 
creasingly well approximated by the form (A^)) ~ CV 2 
[HI [29] f° r some constant C. Therefore, 



dAV(A\V)e 



-jA 



(A3) 



up to a prefactor that can be absorbed in f2. Inserting 



(A3 1 into (All, integrating by parts, and simplifying via 
( ]Al2| yields^]). Though ^ is an approximation, its 
accuracy increases with growing system size. 



In the thermodynamic limit N c /N — 2 — > ^Az, where 
Az = z — Zi so . Therefore, as we must have either a — > 
or 7 — > in the thermodynamic limit, (|5| predicts 
the limiting values a(p) — \Az (7 — * 0) and 7(a) = 
|Az (a — * 0), consistent with the insets of Fig. [7] Note 
that a or 7 vanish, and hence statistics can no longer 
be normalized, when z | z iso . This is because the force 
network ensemble vanishes at isostaticity. 
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